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Abstract 

An ab-initio method for determining the dynamical structure function of an interacting many- 
body quantum system has been devised by combining a generalized integral transform method with 
Quantum Monte Carlo methods. As a first application, the coherent and, separately, the incoherent 
excitation spectrum of bulk atomic 4 He has been computed, both in the low and intermediate 
momentum range. The peculiar form of the kernel in the integral transform of the dynamical 
structure function allows to predict, without using any model, both position and width of the 
collective excitations in the maxon-roton region, as well as the second collective peak. A prediction 
of the dispersion of the single-particle modes described by the incoherent part is also presented. 
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Many important physical properties of matter, from viscosity to magnetic susceptibil- 
ity, are closely related to the underlying microscopic dynamics. The computation of such 
quantities is routinely performed for classical systems. Moreover, thanks to the growing 
computational capabilities, we currently have a number of methods capable of a true ab- 
initio treatment of the ground state of quantum many-body systems. However, the detailed 
study of dynamical properties has been yet elusive. 

In this letter we will focus on the problem of extending in a consistent and reliable 
way a class of Quantum Monte Carlo (QMC) algorithms to make it possible to determine 
the Dynamical Structure Function (DSF) for a generic quantum many-body system and a 
generic excitation operator. Such an extension is based on the use of the Integral Transform 
(IT) technique with generalized kernels. The amount and quality of information that can 
be extracted by the proposed scheme is illustrated in the application to the study of the 
coherent and incoherent density excitation spectrum in 4 He. 

At present, QMC calculations provide benchmark results for the study of a huge variety 
of many-body systems ranging from quantum chemistry, physics of ultra-cold gases, and 
nuclear physics. The most dramatic limitation of QMC methods is their inability to treat 
dynamical properties in a symilarly reliable way. This failure is essentially due to the fact 
that QMC works in imaginary-time rather than in real time. This implies that quantities 
that do not directly translate into imaginary-time language, when analytically continued to 
real time, are affected by statistical noise that can be hardly reduced, making calculations 
unfeasible. 

The mainstream approach to the problem is to attempt a numerical inversion starting 
from the imaginary-time autocorrelation function, i.e. the Laplace transform of the DSF. As 
is well known, such an inversion is an exponentially ill-posed problem [l| and sophisticated 
regularization techniques are needed to correctly extract the physical information. Due 
to the infinite range of the kernel function, the extraction of such informations is extremely 
delicate, and many details of the structure of the excitation spectrum can be missed. As an 
example, one of the most powerful inversion schemes, the Maximum Entropy Method % 5|, 
cannot resolve the (measured) double peaked structure of S(c{, u) in superfluid 4 He 6(, 
corresponding to a higher energy collective roton mode. In a recent paper this structure 
was eventually resolved inverting the imaginary time correlation function in a Path Integral 
Ground State calculation by using a falsification method based on a genetic algorithm 7]. 
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For strongly interacting few-body systems the problem of computing various DSF of den- 
sity and current excitations is solved by using a generalized integral transform approach, i.e. 
the Lorentz Integral Transform (LIT) method |8|, |9|. The success of this approach, applied 
in nuclear physics, lays in the specific choice of the Lorentzian function as kernel of the IT. 
On the one hand this choice allows to calculate the transform with bound state techniques, 
even when the response is defined in the continuum, however avoiding its discretization. 
On the other hand, and most important, the fact that the kernel is a representation of the 
5-function allows for a reliable and stable inversion. So far, the application of this technique 



has been limited to a small number of particles (up to N=6 10j and 7 due to the huge 
computational costs of the diagonalizations needed to calculate the LIT. In the following we 
will discuss how to extend the idea of the LIT to many-body systems, by developing a QMC 
equivalent. 

At zero temperature the contribution to the response of a system of interacting particles 
due to a perturbative probe transferring momentum q and energy uj to it, can be expressed 
using a spectral representation 

Sofau) =EJ(^|0(q)|v&o)| 2 5(^-c) (1) 
= (* |Ot(qM#-w)d(q)|*o) (2) 

where j^o) is the ground state of the system, \^ u ) are the final states of the reaction, O is an 
excitation operator, 5(H — u) is the spectral-density of the hamiltonian and the summation 
is extended to all discrete and continuum spectrum states in the set. 

The cost of a direct calculation of ^(q, u>) becomes rapidly prohibitive as the number 
of particles or the energy transfer u increases. The latter problem is due to the fact that to 
account for continuum states one would need to solve the many-body scattering problem. 

One can instead consider an integral transform of ^(q, u) with a generic kernel K(a, u) 

$(q, a) = J K(a, u) S 6 (q, u) dco. (3) 
The substitution of the expression ([1]) for S^q, u) yields: 

$(q,cx) = ^2^o\dK^u)K(a,co)(^ u \d(q)\^ ) 

V 

= (* \6i(q)K(a,H)d(q)\* Q ) (4) 
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Equation (j3J) can be viewed as a generalized sum rule which depends on a continuous pa- 
rameter a. Provided that the kernel and the excitation operator have suitable analytic 
properties, the right hand side of Eq. (j3J can be calculated using bound-state type methods. 



This is the case both for the Stieltjes kernel 



13, 



18], and for the Lorentz kernel. However, 



while in the former case the inversion of the transform is as problematic as in the case of 
the Laplace kernel, in the latter case even a rather simple regularization procedure allows 
to obtain accurate and stable results. The reason can be easily understood. In the case 
of the Laplace or the Stieltjes kernels the information about ^(q, u) in the u domain is 
spread in a large a domain. On the contrary the Lorentz kernel, as well as any function that 
is a ^-function representation, keeps that information in an arbitrarily narrow a domain, 
governed by the width of the kernel. (In the 5 function limit of the kernel no inversion 
would be needed!) 

The number of 5-function representations is large. However, very few have a practical 
implementation. In the past the use of Gaussian kernels has been investigated in different 
fields from condensed matter ([ijj], [3]) to non perturbative QCD {"hi Tisj] with limited 
results. Here the idea is to recast one possible 5-function representation in the imaginary- 
time propagation language, typical of QMC methods, in a similar way as was suggested in 
Ref. |10|. 

Consider the following family of integral kernels built out of the so-called Sumudu trans- 
form: 

K(a,co) = -\e^ -e-^] P , (5) 
a 



where 



ln\b] — ln\a] ln\b] — ln\a] , 

A* = — 4 -a; " = — h ~b, 6 

— a — a 



and the parameters P, a, b are integer numbers with b > a. The normalization constant iV 
is a function of P, a, b such that J duK(a, u) — 1. 

The efficiency of this transform is displayed by the fact that the kernel function converges 
to a delta function 5(uj — a) in the P — > 00 limit, independent on the choice of a and b. For 
a finite P the kernel aKp(a,u) is still centered around u = a but has a finite width that 
depends on P. These properties make the choice of the resolution and the energy range of 
interest extremely flexible, similarly to the Lorentian kernel. 

Using a binomial expansion, and rewriting powers as exponential functions, leads to a 
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more transparent form of the kernel: 

K P {a iU ) = -Y \\ (-l)V ln ^^ p+ ^. (7) 

By operating the substitution uj — > H according to Eq. (J3J), we are lead to a simple linear 
combination of imaginary-time propagators (h=l), taken at imaginary-time points Tpk = 
bi(b/a)[r^-P + k]/a. Projection QMC methods are all based on the implementation of such 
an imaginary-time propagation. The underlying idea is to solve the corresponding integral 
equation: 

y(R, T ) = J dR'G(R,R',T)^(R',0). (8) 

This can be achieved, for instance, sampling a representation in coordinate space of the 
Green's function G(R, R', r) to propagate a set of configurations representing in turn an 
expansion of the function ^(R,t) in eigenstates of the position operator (Diffusion Monte 
Carlo methods). Alternatively, it is possible to break up the Green's Function in a product 
of short time propagators in coordinate space: 

r) = J dlH... dR n G(R, R n , At) x 

xG(R n , R n ~\ At) ■ ■ ■ G{R", R', Ar)ip(R', 0), 

with t = uAt. This formulation is implemented in the so-called Path-I nteg ral Ground 
State methods [23!], and in the Reptation Monte Carlo (RMC) algorithm [Tjj], where the 
whole path {R, R', R", . . . , R n } is sampled from the product of the short-time propagators 
G, possibly modified with the use of a suitable importance function to be determined 
in a variational calculation. The fact that estimating $(er) reduces to the computation of 
an imaginary time correlation function makes this second formulation more convenient and 
straightforward. In particular, path-based methods yield estimates that never depend on 
the (necessary) importance function used to improve the convergence of the calculation. 

In order to evaluate the transform ([3]) within a QMC approach, we need to compute the 
imaginary-time correlation function, and then construct the corresponding linear combina- 
tions. A small width of the kernel might be achieved either using a large value of P or 
by reducing the value of the ratio b/a. In both cases this would require the computation 
of the imaginary time correlation function for long imaginary time. This might indeed be 
a serious complication. However, a kernel that has a width smaller or comparable to the 
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distance between the typical structures of the DSF is in principle sufficient to extract useful 
informations from the transform. We have also confirmed this in numerical tests. In the 
application described below we used the kernel aK P (a,u) with the typical values of P = 2, 
o = l, and 6 = 2. 

As a first benchmark application to a realistic physical problem, we have considered the 
density excitation response in bulk atomic 4 He at T = 0. The system is modeled as a box 
with periodic boundary conditions, containing N = 64 or N = 125 4 He atoms interacting 
via the HFDHE2 pair- wise potential ( 22j, 21]), which quantitatively reproduces the binding 
energy of bulk 4 He up to the freezing point, effectively including three-body contributions. 
Calculations are performed at the experimental saturation density (n = 0.02186 A -3 ). The 
density excitation operator is defined as: 

N 

6(q) = p = X> iq ' ri , (10) 

8=1 

and the transformed DSF in Eq. ([3]) becomes 

N 

As it is customary in neutron spectroscopy one can distinguish the contribution coming from 
the so-called coherent part, given by the terms with i ^ j, related to collective excitations, 
and an incoherent part with i — j that essentially picks up contributions from single particle 
excitations. We have obtained results for both the full and for its incoherent part, in the most 
studied region of the spectrum: the low momentum phonon-maxon-roton part q ~ 0.3 4- 
3.5 A -1 . Computations have be en p erformed by means of a Reptation Monte Carlo (RMC) 
algorithm, as described in Ref. [19 1. The variational importance function includes two- and 
three-body correlations expanded in a basis set [20] and optimized using a variational Monte 
Carlo procedure. Ground state properties are well reproduced: the ground-state energy per 
particle is e^ MC = —7.23 ± 0.01 K, in good agreement with previous calculations using the 



same potential 19], and with the experimental value e^ v = —7.17 K. The static structure 



factor S(q) is consistent with experimental data and previous calculations ( 19], 20] ). 

Turning to the result on 5*(q, u), the striking difference between the estimate obtained 
inverting the transform with the Laplace kernel or the one in Eq. (jHJ) can be seen in Fig. 
where we compare the results of the inversion obtained from RMC data with both kernels. 
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FIG. 1: A typical result for the response function (Q, ui) in liquid 4 He. Points with errorbars are 
experimental results at Q = 0.4A -1 at T = 1.34K[26j]. The continuous line is the result obtained 
using a Laplace kernel. The dotted line is the result computed by using the Sumudu kernel. 
Theoretical calculations refer to Q = 0.44A -1 , value determined by the size of the simulation box. 

Apart from the small shift of the peak due to the 0.04 A -1 difference in the momentum 
transfer (the momenta are limited by the discretization imposed by the use of a finite sim- 
ulation cell, here L=14.306A ), the new kernel permits to retrieve the information on the 
second peak and gives a much more realistic height and width of the one-phonon peak. 
A few comments are necessary concerning the methods used to invert the transform. 



We have used both the Entropy Maximization Maximum Likelihood (EMML) 25( and the 
Simultaneous Algebraic Reconstruction Technique (SMART) 25]. The bars in the figures 
indicate the computed widths of the excitations. Both the peak position and the linewidth 
are robust with respect to the inversion method used. We found that for q < 2.4A -1 both 
methods converge to the same solution. In Fig. [2] we have plotted the excitation spectrum 



obtained using the new transform. The experimental low-lying part [26( is extremely well 
reproduced up to q ~ 2.6A -1 , where the dispersion does not bend over around 2A but 
continues to raise. In this region, however, the EMML and SMART give different results 
indicating that the statistical uncertainty in the QMC data is too large to allow a consistent 
reconstruction. The two-phonon branch is clearly visible and well resolved. As it happens 
in Ref. it only qualitatively compares to the experiment. However, it should be noted 
that at T = the peak corresponding to the collective excitation should be substantially 
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FIG. 2: Dispersion of the collective modes in liquid 4 He at equilibrium density and T = 0. Points 
with errorbars are the computed values. Errorbars are estimates of the width of the peaks. + and 
x are the corresponding experimental data from Ref. 24(] at T = 1.1K 




0.5 1 1.5 2 2.5 3 3.5 4 

FIG. 3: Dispersion of the first and second peak of the incoherent DSF computed by means of the 
Sumudu integral transform. Empty squares are simulations performed in a simulation box with 
N=64 atoms, stars refer to a box with N=125 atoms. The dashed line is the free-particle excitation 
spectrum. A is the roton gap, and the lines at energy A and 2A are drawn as reference. 



27 
26]. 



. Therefore, the 



narrower. An estimate of the intrinsic peak widht is Au ~ 5 x 10 _4 K 
experimental width is essentially due to the resolution of the apparatus 

The dispersion of the peaks of the collective modes, displayed in Fig. 2, was obtained 
combining two simulations at Ar = 0.002i^ _1 and respectively obtaining a mesh 

separation less than 0.5 untill about 40K. In order to obtain meaningful results in the 
high energy regime, a large collection of RQMC data taken with different imaginary-time 
steps is needed in order to increase the sampling points. Due to this technical difficulty, 
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at present we have not performed an exhaustive research in the high momentum-transfer 
limit. However, some preliminary calculations show that the spectrum has the expected 
approximately free-particle like behavior, and that for Q QA^ 1 and above the incoherent 
part of the Dynamic Structure Factor accounts for the total scattering. 

Indeed, the most useful feature of these calculations is that the resolution is good enough 
to allow for separatly computing the incoherent part of the full response function, in order 
to study single-particle excitations. 

In Fig. 3 we have plotted the calculated excitation spectrum of single-particle excited 
states. The spectrum shows at least two distinct branches. A lower energy excitation 
starts from Q ps and propagates with a velocity resulting slightly higher than the 

superfluid critical velocity ve/vc ~ 1.57. A second branch can be observed starting at an 
energy slightly below two times the roton energy, tending asymptotically to the free particle 
spectrum. Interestingly enough, the lower energy branch crosses the collective excitation 
spectrum exactly at the roton minimum, thereby reinforcing the picture of the roton as a 
single particle excitation of an atom exiting the superfluid. The behavior of these single- 
particle excitations might be significantly affected by the quantum many-body correlations 
induced by the particle-particle interaction. An extensive experimental analysis of the single 
particle excitations properties in superfluid 4 He is unfortunately not available. 

We have proposed a method to extract reliably well resolved spectra from numerical 
calculations implementing imaginay time propagation of an initial state, such as DMC or 
RMC. Computations might be easily extended to the T^O case by using standard PIMC 
methods. The application to the study of the collective and single-particle excitations in 
4 He shows the robustness and the higher resolution power of this technique. The limit to 
the accuracy of the spectra is in principle limited only by the available computer power 
available. 

We thank W. Leidemann, V. Efros, G.V. Chester, M.H. Kalos, J.Carlson and S. Gandolfi 
for useful and stimulating discussions about this subject. Calculations have been partially 
carried out on the supercomputer AURORA in the framework of the AURORA-Science 
project funded by INFN and the Bruno Kessler Foundation (FBK). Authors are members of 
LISC, Interdisciplinary Laboratory for Computational Science, a joint venture of FBK and 
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